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1. Introduction 

Understanding interacting multiparticle systems is a central problem in many areas 
of physics, including condensed matter theory, quantum information processing, and 
complexity science. The difficulty of this problem arises from correlations between 
the particles: Typically, probability distributions of states of two or more particles do 
not factorize, hence the description of interacting systems requires significantly more 
parameters than the description of the same number of non-interacting particles. This 
makes the analysis of multiparticle systems challenging, but it also leads to novel and 
interesting phenomena such as quantum entanglement and classical complex behaviour. 

There are many approaches with which to characterize the complexity and correlations 
of multiparticle systems. In quantum information theory, much research has focused 
on entanglement measures P-[3], but also on other forms of quantum correlations [HIS]. 
Similarly, a number of different measures of complexity have been introduced for classical 
complex systems P4TT]. 

An approach which can be used to measure complexity in the classical as well as in 
the quantum domain makes use of interaction structures [T214T5] . In this approach one 
asks: Given a physical system, can its stationary state or density operator be viewed as 
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a thermal state of an interacting particle system with A;-body interactions only? If this is 
not the case, then how far is the state of the system from the space of all thermal states 
with fc-body interactions? This distance can then be used as a measure of complexity. To 
measure distances of this type it is useful to make use of the relative entropy, and of the 
underlying geometrical structure known in mathematics as information geometry [T6] . 
As already mentioned, an important feature of this approach is the fact that it can be 
used for the quantum and the classical case: In classical multiparticle systems, each 
particle is assumed to be on one of a finite set of states at any one time, and the global 
physical system is described by a probability distribution over the products of the local 
states. In the spirit of Ref. [12] interactions are described by the Hamiltonian function 
on the state space, where a Hamiltonian is a /c-particle Hamiltonian, if it is a sum of 
functions each of which acts on /c-particles only. It is important to stress though that 
the Hamiltonian used in this approach is not necessarily a physical energy function, it 
is mostly a mathematical object with which to characterise the factorisation properties 
of the system's stationary state. Indeed, the approach is applicable to non-equilibrium 
systems as well, for which there may not be any energy function at all. In the quantum 
case, the physical system is described by a density matrix, and the Hamiltonian is an 
operator acting on the corresponding Hilbert space, but again it need not be that of an 
actual physical system. 

A central problem for the practical application of this approach in the quantum 
case is the calculation of the respective distances. For classical systems with binary 
states efficient algorithms are known [TS l fTT t fTS] . For the quantum case, some analytical 
results for states with a high degree of symmetry have been derived [13]. Moreover, 
an algorithmic approach has been proposed in Ref. [19]. Nevertheless the practical 
applicability of these ideas is not clear at present, and more general results and an 
overarching mathematical theory are as yet missing. 

In this paper, we present several results relating to the computation of the distance 
of a given quantum state from the set of thermal states generated by Hamiltonians 
with /c-particle interactions. First, we show how symmetries of the state can be used to 
simplify calculations of this type. Second, we present a new algorithm for the efficient 
computation of such distances. We also discuss how tools from convex optimization can 
be used in order to compute complexity measures for quantum states. 

Our paper is organized as follows: In Section 2, we introduce the relevant notation, 
in particular the formalism of quantum exponential families, and we explain the most 
relevant existing results. In Section 3 we show that for cases in which the quantum 
state being studied carries a certain symmetry the closest thermal state generated by 
/c-particle Hamiltonians has the same symmetry. In Section 4 we present our algorithm 
and discuss its application to specific examples. In Section 5 we discuss how results from 
optimization theory can be used to study this problem. Section 6 finally summarises 
our findings and we present an outlook on future lines of research. 
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2. Exponential families of quantum states 

In this section we introduce the theory of exponential families of interaction spaces for 
the special case of quantum states. We emphasize that most of the results presented here 
have been derived for the classical case by various authors (see, e.g., Refs. [T21IT61I20] ) 
and for the quantum case mostly by D. L. Zhou [T3l[T^. Nevertheless, since these results 
are rather scattered in the literature, we believe that a more comprehensive presentation 
can be useful. 

2.1. Exponential and Block representation 

We consider systems consisting of n two-level systems (qubits) throughout; the 
generalization to higher-dimensional systems is straightforward. For later convenience, 
we first describe two different ways of representing quantum states of such systems. 

The first possible representation is the exponential representation. It uses the fact that 
any quantum state can be considered as a thermal state of some appropriately chosen 
Hamiltonian. More precisely, any n-qubit quantum state of full rank can be written as 

Qc^p{0) = exp( ^ Oai^...,a„(Tai ® " " " ® (Taj (1) 
ai,...,a„ 

where the indices run from to 3 and the cxj are the Pauli matrices with the convention 
ao = 1, cTi = cr^, (72 = CTy and as = a^. In the following, it will be convenient to use a 
multi-index notation, 

Qexp{0) = exp{^e^aa), (2) 

a 

where cJq, = ® ■ ■ ■ ® aa„, and where a = (ai, . . . , The coefficient Oq of the 
identity o"o = 1*^" in the above quantum states is not arbitrary, as it can be determined 
from the normalization condition tr g^xpiO) = 1- Explicitly one has = where 

V;(6>) = ln{tr[exp(^e,a,)]}. (3) 

Any quantum state, g, of full rank can be written in the form of Eq. ([1]), and one 
can view the exponent in the exponential representation as a Hamiltonian of which g 
is the thermal state. In this terminology, the function ip is up to a sign the free energy 
of statistical ensemble [211423] . It is important to note, however, that the Hamiltonian 
does not necessarily correspond to that of an actual physical system. 

An alternative description of the quantum state is given by the affine representation 
or Block representation. In this representation one writes the state as 

a 

where the coefficients are given by r]a = tr(f)afrc"a)- Here, the normalization condition 
is simply given by ?7o = 1. Note that in the affine representation the positivity of the 
density matrix results in additional restrictions on the coefficients rj; these conditions, 
however, cannot normally be formulated straightforwardly [2^123]. 
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We will now briefly discuss the connections between the two representations of 
quantum states. In order to do so consider two states g and g' of full rank in their 
different representations, g = gcxp{G) = f?aff(^) and g' = gexp{G') = gasiv')- Using the 
standard definitioii|| of the relative entropy between two quantum states, g and x, 

D{g\\x)=ti[g\og^{g)]-ti[g\og^{x)], (5) 
as well as the entropy S{g) = — tr(f)log2 g) of a single quantum state, we then have 

H^)D{g\\g') = - \n{2)S[g,s{ri)] - tr{ 1 [l + [ ^ O'^a^ - } 

where the function (l>{ri) = — ln(2)S'[^aff(^)] = — ln(2)S'(f)) is proportional to the entropy 
of g. With the scalar product r} ■ 6' = J2ayto Va^a this result takes the form 

H2)D{g\\g') = <t>iri) + iie')-ri-0'. (7) 
For the special case in which g = g' this reads 

0(77) +^(0)-ry -0 = 0, (8) 

a result which will become important below. At this point it should be noted that the 
expression of Eq. ([8]) shows that ip{0) and (f){r}) are related by a Legendre transformation. 
More specifically, from Eq. ([8]) it follows that rja = dil){0) / dOa and 9^ = d(f){r]) / drja for 
all a ^ 0. 

Similar structures are, of course, well known in statistical mechanics: a 
thermodynamic ensemble in statistical mechanics is defined by the requirement that 
some observables Ai (e.g. the Hamiltonian) have fixed expectation values (e.g. the 
internal energy U). Maximizing the entropy of the statistical distribution of the 
ensemble under these constraints, the thermal state of the ensemble comes out as 
g ~ exp ( — where the coefficients Aj arise as Lagrange multipliers. It is 

then well-known that the Lagrange multipliers Aj are related to the expectation values 
{Ai) by a Legendre transformation (see page 40 in Ref. [21], or [22||23]). 

Let us finally explain a useful theorem for the relative entropies between three states. 
For the pairwise relative entropies of three full-rank states g, g' and g" one has 

D{g\\g") - D{g\\g') - D{g'\\g") = D{g\\g") - D{g\\g') - D{g'\\g") + D{g'\\g') 

ln(2) L 

- 0(V) - ^{0") + ri ■ 0" + (P{rj') + ij{e') -rj' O' 
1 



ln(2) 



iv - rj') ■ {6' - e") (9) 



and thus 



D{g\\g") = D{g\\g') + D{g'\\g") + ^-^(ry - 77') ■ {6' - 6"). (10) 



X Note that the binary logarithm is used in our definition of D{q\\x)- 
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If the scalar product vanishes, this relation is also called the generalized Pythagoras 
theorem |16j . 

2.2. Exponential families and the information projection 

We next define the exponential families of states generated by fc-particle Hamiltonians, 
objects of this type will be the focus of the work presented here. For a given multi-index 
a = (tti, . . . , an) we write for the weight of a, i.e., the number of factors in the 

Pauh operator cTq, = cTq,^ ® ■ ■ ■ ® 0"^^ different from the identity. In other words the 
weight py(Q;) of a multi-index a is the number of nonzero elements Oj. 

For any 1 < A; < n we then can define the so-called exponential family of thermal 
states of /c-party Hamiltonians as 

Qfc = {^1^' = exp( ^ Oa(Ta)}, (11) 
a:W{a)<k 

SO that Qk is the set of all quantum states, for which the exponential representation 
contains only k-hodj interactions in the Hamiltonian. Since the exponential 
representation is unique and the operators form a basis of the operator space, 
this definition is unambiguous. It will also be useful to write Hk for the space of all 
Hamiltonians containing only interaction terms up to weight k. 

The set Qk represents a manifold in the space of all quantum states (see also Fig. 1), 
a direct characterization is not straightforward. Obviously, one has Qk C Qk+i and 
consequently the exponential families define a hierarchy 

Qi c Q2 c ■ • ■ c Qn, (12) 

where Qn is the set of all states with full rank and Qi the set of all product states with full 
rank. For states g with full rank we will then ask what the minimal order of interaction, 
k, is such that g G Qk- For states which are not of full rank, the analogous question is 
whether the state is in the closure Qk of an exponential family. The introduction of the 
closure of exponential families makes the discussions and results that follow applicable 
to quantum states for which some of the eigenvalues vanish. 

For a given state quantum state g one can then construct the distance from the 
exponential family Qk, and in particular the state in Qk which is the closest to g. This 
defines the so-called information projection: 

Definition 1. The information projection gk of a quantum state g is the element of 
the exponential family Qk which is the closest to g with respect to the quantum relative 
entropy, 

gk = aTgmm^,^Q^D{g\\g'), (13) 

where D{g\\x) = tT[g\og2{g)] —tT[g\og2{x)]- The distance to the information projection 
is then considered as a complexity or correlation measure and is given by 

Dk{g)= inf Digy) = DigUk) (14) 




Figure 1. Illustration of the information projection onto a quantum exponential 
family. Shown are the linear family Aikis) of distributions with the same fc-party 
reduced density matrices as g (blue line), the exponential family Qk of thermal states 
of fc-party Hamiltonians (red curve) and the information projection gk of g onto Q^; 
and g' represents an arbitrary state in Qk- See text for further details. 



In order to carry out the computation of this distance it is useful to have different 
characterizations of the information projection gk- To this end, consider a given state g 
and define the set A4k{Q) of states with the same /c-party reduced density matrices as 

Mk{g) = {g' \ g'A = Qa fo^' all A C {l, . . . , n} with \A\ = k}, (15) 

where gA = tT^{ii,...,i„}\A{Q) is the density matrix which is obtained from g by tracing 
out all qubits except those with indices in A. We note that A^fc(^) is a linear subspace 
of the space of all n-qubit density matrices, as opposed to the exponential families Qk- 
Alternatively, one can also write 

Mk[gi.s{r})] = {&ff(V) \Va = Va for all a with W{a) < k}. (16) 

The following Lemma was first proven in Ref. [Hj and presents three equivalent 
constructions of the information projection. It also shows that the state gk in Definition 
1 is unique. 

Lemma 2. The following conditions on a quantum state gk are equivalent [H]: 

(a) The state gk is the information projection in the sense of Definition 1. 

(b) The state gk is the maximizer of the von Neumann entropy in the set Aikio) of all 
states with the same /c-party reduced density matrices as g, 

gk = a.Tgmax S{g'). (17) 

(c) The state gk is the unique element of the intersection of the exponential family Qk 
with the set Aikio) of all states sharing the same fc-party reduced density matrices as 

{Qk} = QknMkig). (18) 
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The situation is illustrated in Fig. [H Instead of computing by minimizing the 
distance from Q^, one can also look at the intersection of with Aikio), or indeed 
maximize the entropy among elements of the linear family Aikig)- 

Proof. Let us start with the case (b), and show that the characterisation (b) is 
equivalent to that of (a). The linear family J^k{Q) is defined as the set of all states g' 
with the same /c-particle reduced density matrices. This condition is equivalent to the 
requirement that tT{g'aa) = Aq, = tT{gaa) for all a with < k. Now, maximizing 

the entropy under the constraint of given expectation values is a well discussed problem 
in statistical mechanics [21]. The following results are known: The state maximizing 
the entropy is of the form gk ~ ^w{^a-w{a)<k ^a^a\ and the maximum is unique. So 
it is clear that the maximization in condition (b) results in a unique state gk in Q^, we 
only have to show that minimizes the relative entropy. 

In order to show that, consider a third state g' in Qk (see also Fig. 1). We apply 
Eq. (fTOj) to the state g = ^?aff(?7), the state gk = ^'aff(^) = Qcxpi^) [defined via Eq. (fT7|) ] 
and the states g' = f?cxp(^') in Qk, resulting in 

D{g\\g') = D{g\\~gk) + D{gk\\g') + ^^(ry - ^) ■ (0 - 6'). (19) 

The terms in the scalar product with iy(a) < k vanish, as we have gk G Aikio) and 
thus rja = fja for these a. The terms with the terms with 11/(a) > k vanish because 
of Oa = — 0- So has D{g\\g') = D{g\\gk) + D{gk\\g'), which implies that gk, as 
defined in Eq. ( ITTI) . is also the unique state minimizing the relative entropy D{g\\g') 
among all g' E Qk- 

Let us now turn to the characterisation (c). The state defined in (b) is obviously in 
the intersection Qk r)M.k{g), so all we have to show is that this intersection consists of 
only a single state. Let us assume the contrary, so that there are two states and g2 in 
Qk n A^fc(^). Then, applying Eq. (ITOl) and the same argument as above one finds that 
= -D(^ill^i) = D{gi\\g2) + D{g2\\gi). Since the relative entropy is positive semidefinite, 
this implies that gi = g2- D 

2.3. Complexity measures: Definitions and Properties 

As already mentioned, a central topic of this paper is the computation of the distance 
Dk{g) as defined in Eq. (|T4l) . Before presenting our results, it is useful to collect some 
of the properties of the distance measure Dk{-). First, note that Dk can increase under 
local transformations, if /c > 2 p3l[T5]. This means that Dk cannot in a naive way be 
viewed as a correlation measure, and so we prefer to call it a complexity measure. 
Next one can define the degree of irreducible k-party interaction as 

Ck{g) = Dk^iig)-Dk{g), k = 2,...,n (20) 

(where D^ = 0). The quantity Ck{g) describes the extent to which the approximation of 
a state, g, improves, if the allowed interactions in a Hamiltonian increase from [k — 1)- 
body interactions to /c-body interactions. By the generalized Pythagoras theorem, the 
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last definition is directly equivalent to 

Ck{g) = DihUk-i), k = 2,...,n-l. (21) 

Furthermore, writing g = ^?aff(^) and Qk = f?aff(^) = Qexp{0), we have, using Eqs. 
and dHD, 

H2)Dk{g) = ln{2)D{g\\gk) = - ln(2)5(f?) + ii6) - r) 
= -\n{2)Sig)+^iO)-'n-0 

= -H2)S{g) + ln{2)S{h). (22) 

This shows that 

Dkig) = Sigk) - Sig), = 1, . . . , n - 1, (23) 

and consequently 

Ck{g) = S{~gk^^) - Sigk), k = 2,...,n-l. (24) 

These different expressions for Dk and Ck can be useful for investigating the performance 
of numerical algorithms to compute the information projection: Having obtained the 
projections gk, one can compute the distances Dk{g) = D{g\\gk) and the interaction 
measures Ck{g)- The latter can be calculated in three different ways, namely via 
Eq. fl20|) . Eq. fl2T|) or Eq. i^^. If gk is not the correct information projection, these 
three expressions will in general give different values. 

Finally, let us discuss the case k = 1 in some more detail. The quantity Di is also 
referred to as the multi-information [20] or the degree of total interaction. It has an 
expansion into a telescopic sum of entropy differences 

n 

Cto,{g)=D^{g) = Y,Ck{g). (25) 

k=2 

This is an orthogonal decomposition in the sense of the generalized Pythagoras theorem. 

The exponential family Qi consists of all product states (with full rank). The 
projection of a state g onto this family is given by the tensor product of the one-party 
reduced density matrices, 

^1 = ^'{1} ® ■ ■ ■ ® g{n} where g{i} = tr{i,...,„}\{i} g. (26) 

For the other projections there is no such explicit formula. Moreover, the family Qi is 
invariant under local filtering transformations of the form 

^ ^ a = [Fi ® F2 ® ... ® F„]g[Fl ® ® ... ® F^] (27) 

where the Fi are arbitrary matrices, since these transformations preserve the product 
structure [26]. This means that the quantity Di cannot increase under these 
transformations either [15]. The exponential families Qk with A; > 2 do not have this 
property. 
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3. Symmetries 

The computation of tlie quantities Dk{g) and Ck{g) is not straightforward. One may 
therefore ask, whether symmetries of the state g can simphfy the optimization procedure. 
For instance, if g is invariant under the permutation of the first two particles, it seems 
natural that the state t E Qk (and the corresponding /c-party Hamiltonian) which 
minimize Dk{g) share the same permutation symmetry. This symmetry assumption 
seems plausible, and the following Lemma presents a rigorous statement: 
Lemma 3. Let ghe a. quantum state which has a symmetry of the form 

g = UgU\ (28) 

where f/ is a unitary matrix that keeps the set of all A;-particle Hamiltonians invariant, 
i.e. it is a transformation so that the transformed operator UHkU'^ is again a /c-particle 
Hamiltonian for any /c-particle operator Hk. Then the state g^ G Qk as well as the 
Hamiltonian Hk minimizing the distance D[g\\T) have the same symmetry as g. This 
means that one can restrict the optimization to states and Hamiltonians which fulfil the 
condition 

T = UrU^ and Hk = UHkU\ (29) 

respectively. 

Proof. Consider a given g with the symmetry and let H G "Hfc be the Hamiltonian 
of the state r G Qk minimizing D[g\\T). Here, H denotes only the terms which are not 
proportional to the identity in the Hamiltonian, i.e. one has r = exp(if) exp[— ?/;(iJ)], 
where exp[— ^/^(if)] is the normalization [see Eq. ([3])]. The idea is to consider H' = 
{H+UHU'^)/2 and the corresponding r' = exp(if') and to prove that Z^(f)||r') < D{g\\T). 
From the conditions on U it then follows that H' G Hk and the uniqueness of the 
information projection implies that H = H' and therefore H = UHW. 

From Eq. ([6]) on sees that D{g\\T) and D{g\\T') differ only in the contributions from 
the normalizations ip{H) [or tplH')]. In fact, D{g\\T') < D{g\\T) is equivalent to 

iIj{H') = ln{tr[exp(i7')]} < ^(H) = In { tr[exp(i7)]}. (30) 

So it suffices to show tT[exp{A + UAU'^)] < tr[exp(2A)] for arbitrary hermitean matrices 
A. Applying the Golden-Thompson inequality tr[exp(A + B)] < tr[exp(y4) exp(5)] (see 
page 261 in Ref. [27]) we have tr[exp(y4+f/74[/"'^)] < tr[exp(y4)f/ exp(yl)t/^] and it remains 
to show that tT{XUXU^) < tr(X^) for X = exp{A). Taking the spectral decomposition 
^ = Y.k^k\(pk){(pk\ this reads Y,,^iCki\k>^i < A^., where = | ((^fc|[/|<^i) p is a doubly 
stochastic matrix, that is, the row sums and column sums of C equal one. Birkhoff's 
Theorem states that any doubly stochastic matrix can always be written as a convex 
combination of permutation matrices, C = "^kPk^k, where the pk form a probability 
distribution (see page 527 in Ref. [28]). So it remains to show that J2k ^k^n{k) < J2k -^fc 
for an arbitrary permutation tt, but this follows directly from the Cauchy Schwartz 
inequality. □ 
This Lemma can be used in various situations to simplify the calculation of the 
information projection: 
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• Permutation symmetry: If g is invariant under permutation of the particles i and 
j, then this is a unitary symmetry as in Lemma 3, with the unitary flip operator 
U = Fij. This unitary operation also fulfils the other conditions of Lemma 3, 
and one can conclude that it suffices to optimize over Hamiltonians with the 
same permutation symmetry. Exploiting this symmetry can reduce the number 
of parameters in numerical algorithms significantly. 

• Graph state symmetry: In the framework of quantum information theory, so-called 
graph states and stabilizer operators have attracted significant attention [29] . They 
are defined as follows: Consider an n-qubit system, and n observables Qi which are 
tensor products of Pauli matrices and which commute pairwise. For example, for 
three qubits one can take gi = ® ciz ® 1, (72 = 1 ® o"^ ® and gs = ax ® o"z ® o'z- 
One can further consider all products of the gi. This is an the Abelian group with 2^ 
elements (since gf = 1.) and this group is called the stabilizer. One can alternatively 
characterize the stabilizer as all tensor products of Pauli matrices, which commute 
with all gi. 

A graph state \G) is defined as an eigenstate of all gi with eigenvalue +1, that is 
gi\G) = \G). For the three-qubit example above, the graph state is given by the 
well-known GHZ state \GHZ) = (|000) + |111))/V2. Allowing also ±1 eigenvalues 
one obtains a basis of 2" graph states \Gk). 

A quantum state diagonal in the basis of the \Gk) is a graph- diagonal state, and 
such states have been intensively studied [SniElj. These states fulfil 



and since the gi = g] are unitary. Lemma 3 can be applied. One can directly see 
that the only possible /c-particle interaction terms in the Hamiltonian which share 
the same symmetry are just all terms from the stabilizer group which are of weight k 
or less. For small k these are typically very few terms, which simplifies calculations 
significantly. Note that this structure was also observed for special graph- diagonal 
states in Ref. [13], but Lemma 3 shows that this holds for all mixed graph- diagonal 
states. 

• U^"^- symmetry: Another family of states where Lemma 3 can be applied are the 
so-called ^/^"'-invariant states. These states fulfil 



for all possible unitary transformations U on a single particle. For two particles 
these states are the Werner states [32], but also for more particles detailed 
characterizations are known [33l|3l]. Lemma 3 shows that when computing Dk{g) 
for these states, the optimal Hamiltonian has the same ?7®"-symmetry. This also 
implies that the optimal Hamiltonian has no single-particle terms, since the only 
single-qubit operator with this symmetry is the identity. It follows that for arbitrary 
[/■^"-invariant states the mult iinformat ion is simply given by 




(31) 



g = U^^giU^) 



(32) 



Di{g) = Do{g) 



n - S{g), 



(33) 
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where Dq{q) denotes the distance to the maximally mixed state. 

4. Iterative computation of the quantum information projection 

In this section we present an algorithm for computing the information projection. Other 
existing algorithms will be discussed in Section 5. 

4-1. Preliminary considerations 

In this section we will describe an efficient numerical algorithm with which to compute 
the projection of a given n-particle quantum state g on an exponential family Qk. 
To this end one needs to construct g^ € Qk such that gk € Mk{g), i.e. such that gk has 
the same fc-particle reduced density matrices as g. The algorithm we put forward is - 
in spirit- similar to the iterative projection algorithm for the classical case proposed in 
Refs. [T71[18]. We start with the fully mixed state g' = 1./ tr(l) as an approximation 
for gk, at each iteration step of the algorithm the current approximation of gk is then 
improved such as to better match the reduced density matrix of g defined by a particular 
subset of particles. The algorithm proceeds by iteratively going through all such subsets 
of at most k particles repeatedly until convergence is reached. 

Specifically, let g be the state whose projection onto Qk we want to calculate, and let 
T = / tr(e^) be the current approximation of pk, where if is a /c-party Hamiltonian. 
The algorithm is initiated from H = H. 

For a fixed £-party observable A with i < k the iteration step consists of adding a 
term eA to the Hamiltonian H, generating an updated approximation r' of g. The 
amplitude of the modification, e, is chosen such that the expectation of A under the 
density matrix r' improves the match with the expectation obtained under g. In detail 
one has the update 

tr(e-f^) ^ ^ tr(e-f^+^^)' ^ ^ 

such that tr(y4r') matches ti^Ag) as closely as possible. In principle e can be obtained 
by solving tr(Ar') = ti^Ag) directly. In practice this is hard to implement though, 
as tr(Ar') depends non-linearly on e. We therefore resort to the following linear 
approximation 



d 

tr [Ar'ie)] = tr(Ar) + ^ ^ tr [Ar (e)] 



+ 0{e'), (35) 



where one has 



d_ 



Ore 



H+eA 



=0 



H+sA 



tr [At'] = tr A^-^-^^ - tr A^-^-j^ trie^+sAy (36) 



tr(9ee 



H+£A\ 



tr(e^+^^)J L tr(e^+^^)- 
We have here used the shorthand notation de = When evaluating the derivative of 
the matrix exponential in this expression one has to take into account that H and A 
generally do not commute. So it is convenient to use the identity |35] 

le^m = f' ds e^Mit)dM{t)_^-sMit)^Mit)^ .37^ 
dt Jn at 
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valid for a general one-parameter family of matrices M{t). Applying this identity to 
substitute for the derivatives in Eq. fl36l) and carrying out a modest amount of algebra 
one obtains 



^tr[Ar'(.)] 



ds ti [Ae^^Ae-^^e^] - [tr(Ar)]'. (38) 

Jo 



s=o tr(e^) 

Next we apply a further approximation to the remaining integral, and replace it by 
the mean of the integrand evaluated at the upper and lower limits s = and s = 1, 
respectively. This gives 

[' dstT{Ae'"Ae-'^e") ^ - [tr(Ae^A) + tT{A^e^)] = tT{A^e^). (39) 
Jo 2 L J 

This in turn leads to 

tr[AT'{e)] ^ tr(Ar) + e{ tr(AV) - [tr(Ar)]2}. (40) 

Admittedly, this approximation can only be justified a posteriori by the performance of 
the algorithm. Setting tr[y4r'] = tr[y4f)] and using the approximation just obtained one 
finds the following solution for e: 

tr(Ag) - tr(Ar) {A),-{A)^ 
^ tT{AH) - [tr(Ar)]2 A2(A) ' ^ ' 

where we have used the notation (A)^ = tT[Ag] (and analogously for (A)^), and where 
A2(A) = (^2)^ - {A)l is the variance. 



4-2. Description of the algorithm 

In the full algorithm with which to compute the information projection of g onto the 
quantum exponential family Qk, one chooses an orthogonal basis Vk in the space of 
/c-party observables (excluding the identity) and updates the approximation r for each 
A e Vfc in turn. For an n-qubit system one can choose the Pauli operators 

Vk = {r^\l< W{a) < k}. (42) 

The complete algorithm is then as follows: 



Problem: Given an n-qubit state g, compute its information projection g^ onto the 

exponential family Q^. 

Algorithm: 

1. Choose an orthonormal basis Vk of the space of /c-party observables, say these 
observables are Ai, A2, . . . , Am (where M will depend on k). For each element 
Ai G Vk compute the expectation value (Ai)^. 

2. Initialize r = 1/2" as the completely mixed state. 

3. (a) Start with i = 1, and update r according to 

e^+-^- _ {A,)^-{A,)^ 

^ tr(e^) ^ ^ tr(e^+-^0 ^ ^ A2(A,) ' 
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(b) Increment itoi+1 and repeat step 3(a). Once the coefficients for all observable 
in Vk have been updated (i.e for z = 1, . . . , M) goto 4. 

4. If the maximum number of iterations has been reached or a convergence criterion 
is met, terminate, otherwise goto 3. 

When implementing the algorithm, it turns out to be useful to introduce an additional 
parameter u which controls the size of the steps in the space of Hamiltonians, 

^ ^ ^ r = -^-^^-^ (43) 

with e as above. Choosing values u < 1 corresponds to what is known as a successive 
under relaxation scheme [36] and it can improve the convergence properties of the 
algorithm. 

We would like to stress that we do not have a proof that the algorithm converges. For 
the classical case, however, it has been shown that a similar algorithm converges [l7l[T8] . 
One could also think of improving our algorithm by using a better approximation to the 
integral. However, the numerical results shown below demonstrate that the algorithm 
as described here works remarkably well. 



4-3. Test of the algorithm 

In order to test the algorithm just described, we consider Dicke states of four and six 
qubits. These states are given by 

\D^) = ^(|0011) + |0101) + |0110) + |1001) + |1010) + |1100)), 

1 

\Dl) = ^=(1000111) + permutations), (44) 



that is, they are a balanced superposition of all terms in the standard basis with n/2 
excitations. Dicke states have intensively been studied in entanglement theory and 
have been observed in several experiments [5]. As a one-parameter family of states, we 
consider Dicke states mixed with white noise. 

Before applying our algorithm, it is useful to discuss the symmetries of these states. 
First, the states g^^\p) are symmetric under the exchange of particles. Consequently, 
the information projection shares the same symmetry. There are, however, additional 
symmetries: If we consider the operators 

G^") = " and G^") = af", (46) 
then it directly follows that = ^^"Hp) for 

a = X, z, which implies 

that the information projection Qk and the corresponding minimising Hamiltonians 
have this symmetry as well. Note that g^'^\p) is also symmetric under the product 
G^^ = G^^G^r\ but this is not an independent symmetry. 
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Figure 2. Complexity measures for Dicke states mixed with white noise (left: the 
fom'-qubit case, right: the six-qubit case). See text for further details. 

These symmetries reduce the number of parameters aheady significantly. For instance, 
there are no single-particle Hamiltonians, which are invariant under these symmetries, 
so we have Di{g) = Dq^q) = n — S{g), as for the [/'^"■-invariant states discussed above. 
Similarly many possible interaction terms can be discarded from the outset for higher- 
order interactions. Applying our algorithm generates the data shown in Fig. |2j We have 
here chosen an under-relaxation parameter of a; = 0.5 for Q^^\p) [w = 0.1 for Q^^\p)], 
although u = 1 yields similar results when convergence is reached. We typically run 
the algorithmic scheme for up to 100 iterations [500 iterations for g^^\p)], or until a 
convergence threshold is met, and we report the minimum distance reached over this 
number of iterations. It must be stressed though that our numerical results are estimates 
of the respective distances, we cannot exclude that the precise quantitative results have a 
remaining dependence on parameters of the algorithm, such as the relaxation parameter 
u, the maximum number of iteration steps, or the precise convergence criterion. 

For the four-qubit case, one finds that the measures D2 and D3 coincide. This can be 
explained as follows: Let us consider the information projection ^2 = exp(if2) and the 
corresponding two-particle Hamiltonian H2 G 7/2- As already mentioned above, H2 does 
not contain any single-particle term. For the two-particle terms there are also not many 
possibilities, in fact, the only possible terms are /i^'^ = cTq, ® cTq, 1 ® 1 for a = x, |/, 2;, 
and permutations thereof. Using the power series of the exponential function, one finds 
that g2 = exp{H2) has no three-body correlations in its Bloch representation, that is, 
tr[(crj (g) CTj ® Cyfc ® 1)^2] = for any choice of i,j,k G In other words, one 

§ The detailed proof is the following: H2 and ^2 obey the symmetry defined by the Ci"-* , which implies 
already that most of the three-body correlations in the Bloch representation vanish. The only terms 
which are not forced to be zero are expectation values of K^^'^ = tTa; (8) cry (8) ct^ (g) 1 or permutations 
thereof. However, if exp(iJ2) is written as a power series, the term K^^^^ does never occur as a product 
of the two-qubit terms /i^-' . The reason is that if the product of the single-qubit observables in K^^") 
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can state that the maximizer of the entropy in the hnear family A^2 has no three-body 
correlations in its Bloch representation. 

On the other hand, let us consider the three-body reduced state of Q^'^^p). This state 
is given by 

Qi23=P^ + - |000)(000| - |111)(111|) (47) 

and one can directly check that this state also has no three-body correlations in its 
Bloch representation. But this means that the maximizer of the the entropy in the 
linear family A^2 has the same reduced three-particle density matrices as Q^^\p), so it 
is also an element of the smaller linear family Ai^. From this D2[0^^\p)] = DsIq^^^p)] 
follows. 

For the six-qubit case, we find that D2[g^^\p)] = D^lg^^^p)] and D/i[g^^\p)] = 
D5[0^'^Kp)] and this can be understood in the same way. In order show that D2 = -D3 
one first finds that the reduced three-particle state of g^^^ (p) does not contain any three 
body correlations, then the argument is the same as for the four-qubit Dicke state. 
For D4 = one has to consider the reduced five-particle density matrices. Due to 
the symmetry, the only relevant correlators are = ax ® Cz o'y ® ay ® ay ® 1, 
K2 = <7x ® <yz ® <yz ® <yz ® <7y ® ^ and = ax ® CFx ® CFx ® CTz ® <7y ® ^- Their mean 
values vanish in the state g^^\p)^ and then the proof can proceed as before. 

Besides these examples, we have tested our algorithm for a variety of four- and five- 
qubit states, and previous results [T3l[T9] we easily reproduced. This shows that the 
algorithm presented above is a useful tool for computing the complexity measure for 
states up to six qubits. 

5. Other algorithms 

In this section, we will first describe another algorithm, which has been proposed to 
compute the complexity measure [TH] . Then, we will discuss how other methods known 
in numerical optimization can be used for this problem. 

5.1. The algorithm of Ref. IW^ 

In Ref. [19] D.L. Zhou has proposed an algorithm for computing the information 
projection and has presented examples up to five qubits. The idea of the algorithm 
is as follows. 

First, one considers the information projection g^ onto the exponential family and 
its logarithm log(^fc), which is effectively the generating fc-particle Hamiltonian. Then, 
according to Lemma 2, g^ obeys the following conditions: (i) First, the mean values of 
Pauli matrices aa in the state gu equal the mean values in the state g, if the weight 

is taken, the result is (ct^) • {<Jy) ■ (cr,) • (1) = il, which is a non-herniitian operator. If an arbitrary 
product of the K^J is considered, and then the product of the single-qubit observables is taken, the 
result is an hermitean operator, since any Pauli matrix occurs an even number of times in the total 
product. 
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W{a) of is smaller or equal k. This is nothing but the condition that the reduced 
/c-particle states of Qk and g are the same, (ii) Second, the mean values of Pauli matrices 
cTq in the Hamiltonian log(^fc) vanish, if the weight M^(q!) of cTq is larger than k. This 
is the condition that the generating Hamiltonian contains no higher-order interactions 
than /c-particle interactions. 

These conditions lead to 4"^ equalities for g^, and from Lemma 2 it follows that there 
is a unique solution to all these equalities, which is the desired Qk- Due to the occurrence 
of the logarithm, however, the equalities are highly nonlinear, and a direct numerical 
solution is not straightforward. Therefore, as explicitly stated in Ref . [19] , one needs an 
initial guess of an initial value of Qk for solving them. For example, in Ref. [19] curves 
like the ones in Fig. 2 have been computed iteratively: Initially, one solves the problem 
if the state is completely noisy {p = 1), then one uses the solution as an initial value 
for solving the nonlinear equations for decreasing p ^ p — e and so on, until p = is 
reached. We stress that no such procedure is needed for our algorithm, then data points 
shown in Fig. [2] are obtained independently for the different values of p. 

5.2. Convex optimization approaches 

The algorithm described in Section 4 searches for an approximation of the information 
projection by an iteration within the exponential family Q^, which is a highly nonlinear 
manifold. In Lemma 2 it was established that the information projection qj. can also 
be characterised by a maximization of the von Neumann entropy S{q') over the linear 
family M.k{Q)i given by the density matrices with the same /c-particle reduces density 
matrices. The advantage of this formulation is that the problem becomes an instance 
of convex optimization, namely the minimization of the convex function —S{g') over 
the convex set Aiki^)- Note that convex optimization problems are well-studied, for an 
overview see Ref. [37] . 

From this structure it is clear that no local minima exist, i.e., if —S{g') < —S{t) for 
all r G Aikio) with — r|| < e where s > 0, then g' attains the global minimum. 
This makes a numerical solution particularly tractable. In the language of convex 
optimization, the problem reads 



For such problems one can construct algorithms where the optimality of the solution is 
guaranteed. 

A problem which is related to the one discussed here was studied by Teo and coworkers 
in the context of quantum state estimation theory [38]: They consider the situation 
where in a quantum experiment the observed frequencies f of measurement outcomes 
(described by positive operators Ei) are sampled from a probability distribution P[g'] 



maximize: t 



subject to: S{g) > t, 



tT[{g — g)aa] = for all a with < k, and 

g' > 0. 
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from an unknown quantum state g' (that is, the probabihties are computed according 
to Pi = tT{Eig')). Among the solutions which maximize the log-hkehhood log{L{g'\i)) 
[which is proportional to — i5(f ||P[^)'])], they propose to use the state with maximum 
entropy. For the solution of this problem they introduce an iterative algorithm that 
eventually approaches 

g* = lim argmin[- log(L(^'|f)) - XS{g)], (49) 

A\0 gl 

where the minimization is performed over all g' > 0. 

In order to see the relation to our problem, consider a probability distribution 
P[r]i = tr(i?jr), where Ei are positive semidefinite operators with "^^Ei = 1 and 
span{i?j} = spanAikig)- This means that knowledge of the probabilities P[t] is 
equivalent to knowledge of the reduced /c-particle density matrices. If we now replace f 
by P[g], then the optimization in Eq. (l49l) is indeed closely related to our optimization 
problem: One can view Eq. ( H9l) as a maximization of the entropy >S'(^') where the log- 
likelihood term [being proportional to i5(P[f)]||P[p'])], serves as a barrier term, forcing 
g and g' to have the same reduced fc-particle states. Such constructions are known from 
interior point methods for convex optimization [37] . 

While the methods presented in this section rely on established numerical algorithms, 
we observed that computing the information projection by solving a convex optimisation 
problem actually requires more resources than the algorithm we propose in Sec. 4. This 
is probably due to the fact, that the algorithm of Sec. 4 exploits the structure of the 
problem in a better way. 

6. Conclusions 

In summary we have used concepts from information geometry to characterise the 
complexity of multiparticle quantum states. Specifically, we considered the distance of a 
given n-particle density matrix from the space of all thermal quantum states generated 
by Hamiltonians with fc-particle interactions. We have shown how symmetries can be 
used to simplify the calculation of the resulting complexity measure. Furthermore, we 
have proposed a new algorithm to compute this measure. This algorithm, we think, is 
computationally more efficient than existing approaches, and in particular we are able 
to compute the above complexity measures for selected six-particle states. 

There are several follow-on problems requiring further attention. First, the complexity 
measure is not yet fully understood, and several interesting open questions remain, for 
example, which are the states with a maximal distance Dk from a given exponential 
family? How is the complexity measure related to known entanglement measures or 
correlation measures? Second, it would be interesting to study this measure in specific 
situations. For instance, for a given n-particle spin model with two-particle interactions 
only one may consider the reduced states of some of the particles and ask, whether 
they still can effectively be described as thermal states of a two-body Hamiltonian, or 
whether higher-order correlations are present. We expect that this may for example be 
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of interest in models of quantum spin chains at or near quantum phase transitions. It is 
known that entanglement measures can be used to study such critical phenomena, and 
so a systematic exploration of the complexity measures we have proposed here cannot 
only help to understand quantum phase transitions better, but also to relate different 
measures of complexity and correlation. 
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